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Abstract 

The finite element solution of two-dimensional anisotropic diffusion problems is considered. 
A Delaunay-type mesh condition is developed for linear finite element approximations to satisfy 
a discrete maximum principle. The condition is shown to be weaker than the existing anisotropic 

1 ^ non-obtuse angle condition. It reduces to the well known Delaunay condition for the special 

case with the identity diffusion matrix. Numerical results are presented to verify the theoretical 
i-^ findings. 
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^ 1 Introduction 

We are concerned with the linear finite element (FEM) solution of the two-dimensional anisotropic 
O diffusion equation 

Tl -V-(DVu) = /, in Q (1) 



*^ subject to the Dirichlet boundary condition 

u = g, on d^l (2) 

where $7 £ is a connected polygonal domain, / and g are given functions, and B = Ili{x,y) is 
the diffusion matrix assumed to be symmetric and strictly positive definite on Q. This boundary 
value problem (BVP) is a model of anisotropic diffusion problems arising in various fields such 
as plasma physics [151 El [13 EH ESI EH] , petroleum reservoir simulation [H |2l [lOl [131 [22] , and 
image processing [6l [T] [211 [Ml [SH] |l3] . A distinct feature of the BVP is that its solution satisfies 
the maximum principle and is monotone when f(x,y) < for all {x,y) G 17. A challenge in the 
numerical solution of the BVP is to design a scheme so that the resulting numerical approximations 
satisfy a discrete maximum principle (DMP). 
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work was supported in part by the National Science Foundation (USA) under Grant DMS-0712935. 
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Development of DMP satisfaction schemes for solving diffusion problems has attracted consid- 
erable interest in the past; e.g., see [iElEliailHllISlEOlESlEZllMlllQlllIllia for isotropic 
diffusion problems where D = a{x,y)I with a{x,y) being a scaler function and [lOl |TT1 [151 ESI 
[IZllSlllMllSSllMllMllMlEQlEIllMlES] for anisotropic diffusion problems where D(x, y) can be 
heterogeneous and anisotropic. For example, Ciarlet and Raviart [9] (also see Brandts et al. [1]) 
show that the linear finite element method for an isotropic diffusion problem satisfies DMP when 
the mesh is simplicial and satisfies the non-obtuse angle condition requiring the dihedral angles of 
mesh elements to be non-obtuse. In two dimensions and for the special case D = /, the condition 
can be replaced by the Delaunay condition, a weaker condition that only requires the sum of any 
pair of angles opposite a common edge to be less than or equal to vr [271 El]- Moreover, Xu and 
Zikatanov [11] show that the non-obtuse angle condition at edges where the diffusion coefficient is 
discontinuous and the Delaunay condition at other places guarantee DMP satisfaction. Recently, 
Li and Huang [28] generalize the non-obtuse angle condition to anisotropic diffusion problems and 
obtain the so-called anisotropic non-obtuse angle condition requiring the dihedral angles of mesh 
elements to be non-obtuse when measured in a metric depending on D. 

The objective of this paper is to extend the Delaunay condition to anisotropic diffusion problems. 
A Delaunay-type mesh condition is developed for the DMP satisfaction of linear finite element 
approximations for those problems. It is shown that the new condition reduces to the Delaunay 
condition for the special case D = / and is weaker than the anisotropic non-obtuse angle condition 
developed in [28]. We attain the new condition by investigating the stiffness matrix as a whole. 
This is different from [28] where only local stiffness matrices on individual elements are considered. 



The main theoretical result is given in Theorem 4.1 



This paper is organized as follows. The linear finite element formulation for BVP ([T]) and ^ 
is given in Section 2. Section 3 is devoted to the description and geometric interpretation of the 
anisotropic non-obtuse angle condition. The Delaunay-type mesh condition is developed in Section 
4, followed by Section 5 with numerical results verifying the theoretical findings. Finally, Section 6 
contains conclusions and comments. 



2 Linear finite element formulation for the model problem 

Consider the linear finite element solution of BVP ([T]) and ^ . Assume that a family of triangular 
meshes {Th} is given for Q. Let 

Ug = {veH\n)\v\dn = 9}. 

Denote by [/^^ C Ug the linear finite element space associated with mesh Th, where is a linear 
approximation to g on the boundary. A linear finite element solution G C/^^ to BVP (jl| and Q 
is defined by 

I {Vv^ fliVu^dxdy= ^ / fv^dxdy, ^ e (3) 

where Uq = U\ with g^ = 0. Generally speaking, the integrals in (3) cannot be carried out 
analytically and numerical quadrature is often necessary. We assume that a quadrature rule has 



2 



been chosen on the reference element K, 



„ mm 

j v{i,r])didr]^\k\^Wkv{hk), ^Wfe = l, (4) 
k=i k=i 

where WkS are the weights and b^s are the quadrature nodes. Many quadrature rules can be used 
for this purpose; e.g., see |12j . An example is Wk = ^ {k = 1,2,3) and the barycentric coordinates 
of the nodes (g, g, |), (g, |, g), and (|, g, g). 

Let Fk be the affine mapping from K to K such that K = Fk{K), and denote = Fxibk), 
= 1, • • • ,m. Upon applying Q to the integrals in ([3| and changing variables, the finite element 
approximation problem becomes seeking € [/^^ such that 

m m 

J2 \K\Y,Wk{'^v'^\KfBib^)Vu%= J2 \K\Y,Wkf{b^)v\b^), yv^eU^ (5) 
K&Th k=l K&Th k=l 

where Vw'^lx and \/u^\k denote the restriction of V^;^ and Vn^ on K, respectively. We have used 
the fact that Vv^\k and Vu^\k are constant in deriving nbh. Let 



nK = ^wkn{bk)- (6) 
fc=i 

Obviously, D/^ is an average of D on K. Eq. ^ can be written into 

m 

^ \K\ {Vv%f Bk Vu^Ik = Yl \K\Y,Wkf{b^) v\b^), yv^ G U^. (7) 
Ken KeTh k=i 

We now express ([T]) in a matrix form. Denote the numbers of the elements, vertices, and interior 
vertices of mesh Th hy N , N^, and N^, respectively. Assume that the vertices are ordered in such 
a way that the first vertices are the interior vertices. Then Uq and can be expressed as 

Uq = span{(/)i, • • • , 4>N,,}, (8) 
i=i j=N^i+i 

where (pj is the linear basis function associated with the j-th vertex, aj. The boundary condition 
([2]) is approximated by 

Uj=g{aj), j = N^i + l,...,Ny. (10) 

Substituting ^ into and taking = (pi {i = 1, N^) in ([t]) and combining the resulting equations 
with (10), we obtain the linear algebraic system 

Au = f, (11) 

where 

A=^ , (12) 

/ ' ^ ^ 
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and / in (12) is the identity matrix of size (N^ — N^)- The entries of the stiffness matrix A and 

(13) 



the right-hand-side vector / are given by 



'i\K 



V0 



>j\K, 



i = l,...,N^i, j = l,...,Ny 



K&Th k=l 



(14) 



The expression (13) can be simphfied. Let u!i be the patch of the elements sharing vertex aj. 



Noticing that V<f>i = for {x, y) ^ Wj, we have, for i ^ j , i = 1, N^, j = 1, N^^ 



\K\ {V(t>i\Kf y<Pj\K + \K'\ {V(t)i\K'7 ^K' y^jW'. 



(15) 



In (15), K and K' denote the two elements sharing the common edge (e^j) connecting vertices 
= af = af and aj = = a^' ; see Fig. 1 




af (af) 

Figure 1: Elements K and K' share the common edge (ejj) connecting vertices af^ (af ) and 
(af)- The angles opposite the edge are denoted by and a^', respectively. The Delaunay 
condition is a.^ -|- a/^ < tt. 



3 The anisotropic non-obtuse angle condition 

In this section, we study mesh conditions under which the linear finite element scheme ([7]) satisfies 
DMP. 

To start with, we introduce some notation. Denote the vertices of an element K by af , a^, a!^ . 
The edge matrix of K is defined as 

EK = [af -af,af -af]. 
Since K is simplicial, Ek is nonsingular [37|. A set of q- vectors (cf. Fig. [2]) can then be defined as 

W^\qf]=E],^, qf = -qf-qf. (16) 
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a 



K 



9f 



a 



K 



Figure 2: A sketch of the q vectors and other geometric quantities for an arbitrary element K. 



By definition, qf is the inward normal to the edge opposite to vertex af (i.e., the edge not having 



O'j^ as a vertex). This orthogonality implies that the (dihedral) angle, a^, opposite to edge Cij can 



be calculated in terms of qf and qj^ as 



a. 



K 



TT — arccos 



rrKW \\^K\ 



(17) 



Moreover, it is known [3l [23] that 

From this relation, it is not difficult to show 
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(18) 



(19) 



where hf is the height of K in the direction of qf or the shortest distance from af to the edge 
opposite to af\ see Fig. [2| 

Now, we are ready to describe the anisotropic non-obtuse angle condition. 



Lemma 3.1 If the mesh satisfies the anisotropic non-obtuse angle condition 



■.K\Ty 



qf<o, y i^j, i,j = 1,2,3, y K e% 



(20) 



then the linear finite element scheme ^ for solving BVP ^ and ^ satisfies DMP. 

This lemma was proven in |28j in any spatial dimension by showing that the stiffness matrix A 
in (11) is an M-matrix and has non-negative row sums. A key step of the proof is to show aij < 
for all i 7^ j, which can be seen to hold from (15), (18), and (20). 

For the isotropic diffusion case, the condition (20) reduces to 

iqffqf<0, y i^ j,i,j = 1,2,3, yKeTh. (21) 



Thus, (20) is a generalization of (21) for a general diffusion matrix. Notice that (21) implies that 
the second angle on the right-hand side of (17) is between tt/2 and vr. Consequently, (21) is exactly 
the non-obtuse angle condition [9j, implying < '/r/2. 

The condition ( 20 ) can be more directly interpreted as requiring the angles of elements to be 
non-obtuse when measured in a metric depending on D. To see this, we first notice that, according 
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to (20), the angle between qf and q^ should be measured in the metric Dj^. Indeed, it has the 



expression 

arccos 

where the Di^'-norm is defined by 



{qffOKqf 

kf IIbk Ikf lb 



V-yGM^ (22) 



Since 



/ {qffBKqf \ / (Digf )^(Diqf ) 
arccos ^ — = arccos j ^ — 

Uk mj WokJ \||D|qf II • iio^qf I 

1 1 

the angle can also be regarded as the one between vectors ISj^qf and D^qj^ in the Euclidean norm. 
Denote the third vertex of K by a^. By definition, qf and qj^ are orthogonal to edges (aj^ — a^) 
and {af — ), respectively; i.e., 



It follows that 



(af-af) = 0, iqff{af-af) = 



Iqff ( D7(af - af ) ) = 0, {Blqff (©^(af - af ) ) = 0, 



indicating that Bl^qf and O^qj^ are orthogonal to {a^ — af) and B^^^ {af — af), respectively. 
Thus, the angle between edges (aj^ — af) and (af — af) in the B]^ -norm and that between qf 
and qj^ in the B^ norm are related by 

arccos — xn + arccos ( ^^^^ ^^^^ ) = n. (23) 




Since (20) means the first angle on the left-hand side of the above equation is between tt/2 and 



vr, we conclude that condition (20) is equivalent to the requirement that the angles of elements be 



non-obtuse when measured in the B>j^ norm. 



It should be emphasized that condition (20) has been obtained by considering only local stiffness 



matrices on individual elements. For the current 2D situation, this means that each term in (15) 
has been required to be non-positive. Clearly, this is too strong since we only need Oij < for i ^ j 
for A to be an M-matrix. For the special case B = /, the Delaunay condition requiring the sum of 
any pair of angles opposite a common edge to be less than or equal to vr (cf. Fig. [T]) is sufficient 
to guarantee < for i ^ j. It is then natural to ask if condition (20) can be weakened and a 



Delaunay-type condition exists for the general diffusion matrix B. This issue is studied in the next 
section. 



4 A Delaunay-type mesh condition 

In this section, we develop a Delaunay-type mesh condition under which the linear finite element 
scheme ([7|) satisfies DMP. The main result is given in Theorem 4.1 Its proof is broken into a series 
of Lemmas. 
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Lemma 4.1 For any element K, 



\K\{Vct>i\K)'^V(t^j\K = -^coi{af^), i^j,ij = l,2,3 (24) 



where afj is the angle between edges eu and ekj , with being the third vertex. 

Proof. This result has been obtained in |14j . For completeness, we give a short proof here. 
Without loss of generality, we consider the case with z = 1, j = 2, and k = 3 (cf. Fig. [2]). From 
(17]), ([18]), and we have 

\K\{V(l)i\KfV(l)2\K = \K\{qffq^ 

= |i^| llqfll • ||qf||cos(^-Qf2) 

cos(ai2)- 



hfhf 

From Fig. [2] it is easy to see 



1 h^h^ 

1^1 - 7^"'2 ll^l - II - ^ . f K\ 

I 2sm(af2) 



Combining the above results, we obtain inequality (24). 



The angle can be calculated in terms of the q vectors as in (17) or in terms of the edge 
vectors as ^ ^ ^ ^ ^ 

af^ = arccos n V t,, u'k ■ (^5) 

\ ll"j "A: II ll"j "fc II / 

The above formula is more desirable if linear coordinate transformations are involved. This is 
because, under a linear coordinate transformation, the edge vectors of K will remain to be the edge 
vectors of the transformed element but in general the q vectors will not. The latter is due to the 
fact that orthogonality between vectors is not preserved by linear coordinate transformations. 

Lemma 4.2 For any element K, 

\K\{VcPi\KfBKVcPj\K = -^^^^|5^cot(a^_^_0, i ^ j, i, j = 1,2,3 (26) 

where a^^_i is the angle between edges eki and ekj (with being the third vertex) measured in 
the metric ^J^^ , i-e., 

= ai'^^o^ lla^-a^ll illa^-a^ll i ' ^ ^ 

Proof. Consider a linear mapping G . K ^ K defined as 

M Ud'M V(x,y)Gi^ (28) 



7 



where K = G{K) and (x, y) and r/) are the coordinates in K and K, respectively. Let = G{af ) 
{i = 1,2,3), Cij = G{eij) {i / j), and V = {{d/d^), {d/dr]))^. Denote the angles of K by aij. It 
is easy to show that e^j's form the edges of K and r]) = t])) {i = 1, 2, 3) form the 

linear basis functions on K. Moreover, 



V 



"K 



1 



Since V(/>j and ^<pj are constant on we have 



^V(^jdet(D|,)(i^dry 
v'det(D^) |i^| (V0,;|^,)^V(^j|^. 



K 



Applying Lemma 4.1 to the last term in the above equation on element K, we have 



\K\{V^i\KfOKV^,\K 



Vdet(D^) 



cotld 



(29) 



From 



and similar formulas for {af — a^), diij can be expressed as 



a 



k ); ll^i ~ 0,k\ 



a 



K K 
fJ-i. 



ai 



arccos 



arccos 



— ttfcll \\dj — dk 



a 



Combining this result with ( 29 ) gives ( 26 ) . 



Lemma 4.3 The entry Uij of the stiffness matrix A, (15), can he expressed as 

A/det(Di^) 



COt(Q^„_/ 



\/det(DKO K' ^ 

cot a „_i , 

2 ^ ^o^^K' 



Proof. This lemma is a consequence of combination of (15) and Lemma 4.2 
Theorem 4.1 If the triangular mesh satisfies 



/ det(Dj^) ^ , 
det(D;^0^ ^'^'"^'^ 



/ det(DxO /f' . 
arccot \ — — -- — rcot a..„_i, 
I V det Di^ ^ 



(30) 
□ 



< vr, for all interior edges aj (31) 



where K and K' are the elements sharing eij, then the linear finite element scheme ^ satisfies 
DMP. 
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Proof. We first show that if the mesh satisfies 



a^.^-i + arccot 



/det(DA'' 



cot(a'^ i) ) < vr, for all interior edges e 



det(I 



(32) 



then the conclusion holds. Indeed, notice that the inequality 

V^det(Dx 



+ COt(Q 



> 



can be written as 



K 



«"„ _i < arccot 



/det(D^ . 
det(DA) 



vr — arccot 



/det(I%0 K' ^ 



which is exactly (32). Then, from Lemma 4.3 we have < for i = 1, ...,Nyi and j = 1, ...,Ny 
if (32) is satisfied. The result also means Ojj < for all i ^ j due to the special structure (12) 
of the stiffness matrix. Following the proof of Theorem 2.1 of [28] we can then show that A is an 
M-matrix and has non-negative row sums, which implies that the linear finite element scheme ([7|) 
satisfies DMP (cf. Stoyan [40j or Lemma 1.2 of [28]). 



Next, it is easy to show that (32) is equivalent to 



K 

a. . „_i + arccot 



I detjBK) K ^ 
det(Dx') ' 



< vr. 



(33) 



As a result, (31) and (32) are mathematically equivalent. 



We now study the mesh condition (31). We first consider the case with constant D. For this 



case, 



Then (31) reduces to 



"ii B-i + "ii B-i - ^' interior edges e 



D, det(Dif) = det(Di^/) = det(D). 



For the special case with D = /, (34) reduces to 



afj + afj < vr, for all interior edges e. 



(34) 



(35) 



which is exactly the Delaunay condition (cf. Fig. [T]). Thus, the mesh condition (31) reduces to the 
Delaunay condition for the special case D = / and is a generalization of the Delaunay condition for 
a general D. 

Next, we consider the mesh condition 

vr 



i^ j, ij = 1,2,3, ^KeTh 



(36) 



for a general matrix-valued function 



!)(x,y). From (23) and (27) it is not difficult to see that 



this mesh condition is equivalent to the anisotropic non-obtuse angle condition (20). Moreover, 
under (36) we have 



a^__i < — , cot(a^^_i) > 0, arccot 



det(I 



'K) 



det(I 



'K' 



cot (a 



K 



< 



vr 
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and similar results for a^'p-i and thus (31) is true. Therefore, (36), or equivalently (20), implies 



(31). In other words, the mesh condition (31) is weaker thantne anisotropic non-obtuse angle 



condition (20). 

Finally, we consider some special cases for (31). It is obvious that (31) reduces to (34) when 



det(DA-) = det(Di^/). In Fig. 
plotted for two cases where t 



the region of (a^^_i , at . „_i 



satisfying the mesh condition (31) is 



le ratio det(D7i-')/det(Dx) is either large or small. From the ngure, 
~ K' 



one can see that when the ratio is large (Fig. 3fa)), a'\^.' ^ should essentially be non-obtuse whereas 

|_f 

a^jj_i can basically be any angle between and vr. On the other hand, when the ratio is small 
(Fig. 3 'b)), the roles of _i and i switch. This observation is consistent with that made by 



Xu arid Zikatanov [H] that the non-obtuse angle condition should be imposed at edges where the 
diffusion coefficient is discontinuous (and thus the ratio det(Dx')/det(Dii') can be large or small) to 
guarantee DMP satisfaction. It is also interesting to observe from Fig. [3] that the DMP satisfaction 
region overlaps with + a^.' > vr. 



(a): v^det(Di^')/det(Dii-) = 100 



(b): v^det(DA")/det(I 



>K) 



0.01 



1.5 2 2.5 3 

alpha(ij, Kprime) 



1.5 2 3.5 3 

alpha(ij, Kprime) 



Figure 3: The x and y axes are a^' i and i, respectively. The DMP satisfaction region 



(satisfying the mesh condition (31)) is below the plotted curve. 



5 Numerical results 

In this section, we present some numerical results obtained for BVP ([T]) and ^ with 

f = 0, g{x,0)=gil6,y) = 0, 

/ 0.5y, for0<y<2 Jl, for < x < 14 

q O,w =< „ and q(xAo) = < 

1 1, forl4<y<16 ' ^ ^ 8 - 0.5x, for 14 < x < 16. 



The diffusion matrix is defined as 

B{x,y) 



500.5 499.5 
499.5 500.5 



This example has a constant but anisotropic diffusion matrix D and a continuous boundary condi- 
tion. It satisfies the maximum principle and its solution stays between and 1. 
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The computation is done with four types of triangular meshes shown in Fig. |4j Meshes (a) and 
(b) are obtained by dividing a rectangle into two triangles using the northwest diagonal line and 
the northeast line, respectively, Mesh (c) obtained by dividing a rectangle into four triangles with 
the intersection toward the northeast corner, and Mesh (d) is a Delaunay mesh (which satisfies the 
Delaunay condition). As mentioned in the previous section, mesh condition (31) reduces to (34) for 
the current example (with constant D). Note that Meshes (a) and (d) do not satisfy (34) whereas 



Meshes (b) and (c) do (cf. Fig. [4]). Especially, Mesh (c) has obtuse elements (with angles greater 
than 7r/2 in the D^-'^-norm). 

Fig. [5] shows the contours of the linear finite element solutions obtained for meshes finer than 
those shown in Fig. |4] One can see that finite element solutions for both Meshes (b) and (c) 
stay between and 1 and show no undershoots and overshoots. This is consistent with Theorem 
On the other hand, both Meshes (a) and (d) lead to undershoots and overshoots in the 



4.1 



computational solutions. Fig. [Hjshows these undershoots and overshoots as functions of the number 
of mesh elements. As the mesh is refined, the undershoots and overshoots decrease very slowly and 
eventually reach a rate 0{N~^'^), where A'^ is the number of elements. 




10 12 14 16 




Figure 4: Numerical example in Section [5j Meshes used in computation. The maximum values for 



a 1 and (q 1 +a 1), respectively, are O.OSvr and l.OGvr for Mesh (a), 0.497r and O.OSvr for 
Mesh (b), 0.5l7r and vr for Mesh (c), and O.DSvr and 1.967r for Mesh (d). 
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(a) 





/ /// ■v// 

























Figure 5: Numerical example in Section [5] Contours of linear finite element solutions. 



6 Conclusions and comments 



In the previous sections we have developed a Delaunay-type mesh condition ( |31[ ) under which the 
linear finite element scheme ([T]) for solving the anisotropic diffusion problem ([T]) and ([2]) satisfies 
DMP. This condition is weaker than the anisotropic non-obtuse angle condition ( |20[ ) developed in 
|28j . It reduces to (34) when the diffusion matrix D is constant and especially to the Delaunay 
condition when D = /. The main theoretical result is given in Theorem 4.1 and verified by numerical 
results. 

It is well known that the Delaunay condition can be satisfied by a Delaunay mesh which can 
be generated through edge swapping from an existing triangular mesh. Moreover, Mlacnik and 
Durlofsky [32] have demonstrated that a properly designed edge swapping procedure can improve 
the monotonicity of finite volume approximations for anisotropic diffusion problems. Clearly, the 



mesh condition (31) can serve as a criterion for designing such a procedure. The development 



of an edge swapping procedure based on (31), the convergence study of edge swapping, and the 



generation of a mesh satisfying (31) through edge swapping may deserve future investigation. 
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(a) : For the type of mesh in Fig. |4]j^a) . (b) : For the type of mesh in Fig. |^ 





10000 
N (number of elements! 



Figure 6: Numerical example in Section [5] Overshoots and undershoots as functions of the number 
of mesh elements. 
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